N2O AF calcs v1, with polarization dependence¶
18/05/24
Version with polarization dependence.
22/03/24
Alignment frame calcs preliminary:
- Load and plot alignment data (ADMs) for all temperatures.
- Compute AF results for given temp (with t downsampling), note default case assumes linear pol, $z$-axis alignment.
- Plot XS and betas, interactive heatmaps and line plots.
Note all plots are interactive, mouse-over for data points, and use toolbars for zoom/selection, and widgets for browsing data types (where applicable).
For general methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
In [1]:
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
Load alignment data¶
Now with ADM class.
In [2]:
from epsproc.classes.alignment import ADM
from pathlib import Path
dataPath = Path('~/ePS/N2O/alignment/Wigner D polyatomic normalized')
ADMtype = '.txt'
ADMin = ADM(fileBase=dataPath.expanduser(), ext = ADMtype)
ADMin.loadData(keyType='setADMsT', normType='wignerDpoly')
In [3]:
# import epsproc as ep
# ep.setADMs
In [4]:
ADMin.plot(width=1000)
# ADMin.plot(xlim=(30,50), width=800) # Auto-stack plus pass args & zoom in on specific region (note slider will reset region with interactive zoom)
Load matrix elements¶
In [5]:
# Plotters
# from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [6]:
# # Scan for subdirs, based on existing routine in getFiles()
# fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
fileBase = Path('~/fock-mount/globalhome/eps/N2O/N2O_valence').expanduser() # Data dir on Jake
In [7]:
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0)
data.scanFiles()
# data.jobsSummary()
Compute AF¶
Define field(s)¶
Using EfiledPol class, set various ellipticities and fields $(E_x,E_y)$. For MF calcs, compute for all polarization geometries corresponding to fields aligned along $(x,y,z)$ molecular axes.
Note default case assumes $z$-axis alignment, and $E_y$ for the field plots here is reoriented along the alignment axis in the AF calcs. See E-field pol docs for details and control over geometries.
18/05/24: debugging for alignment... in this case RX is not used in calcs, so need to ROTATE FIELD TERMS prior to PAD calcs.
In [8]:
from epsproc.efield.epol import EfieldPol
import numpy as np
import matplotlib.pyplot as plt
# Mulitple ellipticities and rotations
# For fields, set with format = [amplitude,azimuth (0,pi), ellipticity (-pi/4,pi/4)]
states = 3
# angles = np.linspace(0, np.pi, states) # For range of angles
angles = np.zeros(states) # no rotation
maxEll = 0.5 # Set max ellipticity to use (0 - 1)
ellipticities = np.linspace(0,maxEll*np.pi/4, states) # Set for 0 - max ellipticity
ell = np.c_[np.ones(states),angles,ellipticities]
labels = (ellipticities/(np.pi/4) * 100).round(2) # Set labels at %age ellipticity
# Set fields
Eell_multi = EfieldPol(ell=ell, labels = labels)
# Plot fields
Eell_multi.plot(figsize=(6,6), draw_arrow=False)
plt.legend(labels, loc='upper right')
Out[8]:
In [9]:
# Optional - set field orientation for PAD calcs. If not set defaults to Ey>z-axis rotation.
# Eell_multi.setOrientation
Eell_multi.setOrientation()
In [10]:
import epsproc as ep
ep.setPolGeoms
Out[10]:
In [11]:
from epsproc.geomFunc.geomUtils import EfieldPolConfig
ep, p, RX, EPRX, EfieldPol = EfieldPolConfig(Eell_multi)
In [12]:
ep
Out[12]:
In [13]:
p
Out[13]:
In [14]:
RX
Out[14]:
In [15]:
EPRX
Out[15]:
In [16]:
Eell_multi.YLM
Out[16]:
In [26]:
%matplotlib inline
# Note some build envs may need mpl inline reset here!
# Plot light fields - input
# Eell_multi.plotSph(dataType='rot', facetDim=['t','Euler'])
from epsproc.sphPlot import sphFromBLMPlot
Itp, fig = sphFromBLMPlot(Eell_multi.YLMrot.sel(t=0).squeeze(drop=True), plotFlag = True,
backend='mpl', facetDim='Euler')
In [ ]:
Itp, fig = sphFromBLMPlot(Eell_multi.YLMrot.sel(t=0).squeeze(drop=True), plotFlag = True,
backend='mpl', facetDim='Euler')
In [20]:
Eell_multi.YLMrot
Out[20]:
In [27]:
Eell_multi.EPRX
Out[27]:
In [28]:
from epsproc.sphCalc import setPolGeoms, setBLMs, TKQarrayRotX
TKQarrayRotX(EPRX, RX)
Out[28]:
In [29]:
EPRX
Out[29]:
In [ ]:
break
TO DO May 2024: E-FIELD ROTATION... calcs below seem to be for xpol case otherwise (XS show opposite alignment dependence to previous calcs)
In [31]:
Eell_multi.YLM
Out[31]:
In [32]:
Eell_multi.YLMrot
Out[32]:
In [43]:
Eell_multi.YLMrot.swap_dims({"Euler":"Labels"}).unstack('BLM').sel(Labels='z')
Out[43]:
In [54]:
# Test EPRX setting from rotated case - currently missing...
from copy import deepcopy
Etest = deepcopy(Eell_multi)
# Test replacing epDict with rotation (subselection)
subselection = Eell_multi.YLMrot.swap_dims({"Euler":"Labels"}).unstack('BLM').sel(Labels='z').squeeze()
Etest.epDict = {'ep':subselection.values,
'p':subselection.m.values,
'RX':Etest.RX}
Etest.epDict
Out[54]:
In [55]:
Etest.calcEPR()
In [56]:
Etest.epDict
Out[56]:
In [52]:
hasattr(Etest,'YLMrot')
Out[52]:
In [51]:
Etest.EPRX.unstack()
Out[51]:
In [30]:
Eell_multi.epDict
Out[30]:
AF calc¶
In [57]:
# Set ADMs to use for calc - set for specific T and downsample
# Subselection, note temperature is set by dataKey here
temp = '10K'
trange = [38,42]
# ADMin.subsetADMs(dataKey = temp, trange=[35,45],tStep = 5)
ADMin.subsetADMs(dataKey = temp, trange=trange,tStep = 5)
# Plot subselection
ADMin.plot(keys='ADM')
# Set to main data structure for calcs
data.data['ADM'] = ADMin.data['ADM']
tpoints= data.data['ADM']['ADM'].t.size
In [58]:
# Basic routine (no grouping)
thres = 1e-5
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-2) # OK for 41 t-points (x2 downsample), but missing some orbs and t-points in output!
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres) # OK for 41 t-points (x2 downsample), still some missing regions for thres=1e-2
# mostly contiguous for thres = 1e-3,1e-4
# EfieldInput = Eell_multi
EfieldInput = Etest # Testing frame rotation...
# With pol...
data.AFBLM(AKQS = data.data['ADM']['ADM'].real, EfieldPol = EfieldInput,
selDims = {'Type':'L'}, thres=thres)
# For thres=1e-3, 3 pol states, 41 t-points, requires ~ 50Gb RAM/10 mins. Some dicontinuities
# May need grouping for larger calcs.
# Thres = 1e-4, 3 pol states, 16 t-points, similar requirements
# Testing 22/05/24 - implementation of field frame rotation for EPRX
# Save full class with Pickle
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
import pickle
outfile = f'N2O_AFBLM_polDep_{temp}_t{trange[0]}-{trange[1]}_tp{tpoints}_{thres}_frameRot_220524.pickle'
with open(outfile, 'wb') as handle:
# data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data.data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data.data with pickle, file: {outfile}')
In [ ]:
# # Load data
# thres = 1e-5
# # infile = f'N2O_AFBLM_{temp}_{thres}_220323.pickle'
# # infile = 'N2O_AFBLM_10K_0.0001_220323.pickle'
# infile = 'N2O_AFBLM_10K_1e-05_220323.pickle'
# import pickle
# with open(infile, 'rb') as handle:
# dataIn = pickle.load(handle)
# data.data = dataIn
Plots¶
NOTE: interactive plots. Roll over plots for values and tool options. For heatmaps use the "box select" tool on the histogram at the side to change the colour mapping scale.
Remarks:
- Some unphysical spikes in values <5eV, these are regions with small cross-sections, and are currently filtered out in these plots.
- General patterns seem to be simlar to R-matrix results, as expected, although absolute Blm values seem larger...? Might be a normalisation difference here, TBC.
Cross-sections¶
In [ ]:
# v2 with pkg code...
# hvTest = data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = False, addADMs = False, XS=True, renderPlot=False, returnPlot=True)
data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = False, XS=True, clim=(-1,2))
In [59]:
# Plot some results
# Example with post-plotter subselection for faster plotting...
# Create HV object, but don't plot
data.BLMplot(dataType='AFBLM', xDim='pol', thres = 1e-2, backend='hv', width=700, ylim=(-1.5, 2.5),
renderPlot=False)
# Replot from Holoviews dataset with selectors
# Note this may need hv or hvPlotters to be loaded
from epsproc.plot import hvPlotters
xDim = 't'
# data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).select(Orb='orb9_P').to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm','pol'])
# data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm'])
data.plots['BLMplot']['hvDS'].to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm'])
Out[59]:
In [ ]:
data.data['ADM']['ADM']
In [ ]:
# from epsproc.plot import hvPlotters
# hvPlotters.setPlotDefaults(fSize = [700,300], imgSize = 600)
data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = True, XS=False,
clim=(-1,2), cmap='vlag')
In [ ]:
data.plots['BLMplot']['XR'].attrs['harmonics']['dimList']
In [ ]:
def testKWARGS(**kwargs):
if 'clim' in kwargs.keys():
print(f"clim = {kwargs['clim']}")
else:
print("Nope")
testKWARGS(nothing='thing')
testKWARGS(clim='thing')
testKWARGS(clim=(-1,2))
In [ ]:
plotObj = data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = True, XS=False,
clim=(-1,2), cmap='vlag', returnPlot = True, renderPlot = False)
plotObj.hist(bin_range=(-2, 2))
# bin_range=(-5, 5)
In [ ]:
break
# hvTest.hist()
# hvTest.overlay(data.plots['BLMplot']['XR'].attrs['harmonics']['dimList']) #.opts(**kwargs)
In [ ]:
abs
In [ ]:
subset = data.plots['BLMplot']['XR'].copy()
# selDims =
import epsproc as ep
ep.matEleSelector(subset.unstack(), thres=thres, inds = {'l':slice(2,20)}, drop=True) #, dims = contiguousDims, sq = sqSelector)
# ep.clea
subset.where(subset.l>0, drop=True)
In [60]:
from epsproc.plot import hvPlotters
import numpy as np
# Set data - use existing plot routine but skip line outs (quite slow, replotted at end of notebook)
data.BLMplot(backend='hv', xDim = 't', renderPlot=False)
# Heatmaps for cross-sections
xrAF = data.plots['BLMplot']['XR'].copy()
Etype = 'Eke'
# Set opts to match sizes - should be able to link plots or set gridspace to handle this?
hvPlotters.setPlotDefaults(fSize = [750,300], imgSize = 600)
cmap = 'coolwarm' # See http://holoviews.org/user_guide/Colormaps.html
# Get cross-sections
XS = np.abs(xrAF.XSrescaled) # Should be in MB
XS.name = 'XS/MB'
hvXS = hvPlotters.hv.Dataset(XS.fillna(0))
# Plot heatmap + ADMs
(hvXS.to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1)
Out[60]:
AF-$\beta_{L,M}$¶
In [ ]:
# Heatmap for l=2,4,6
# NOTE - clipped cmap to physical values with clim=(-1,2), but some spikes.
# BLM = np.abs(xrAF.XSrescaled) # Should be in MB
# XS.name = 'XS/MB'
# Filter out small XS regions?
xrAFfiltered = xrAF.where(np.abs(xrAF.XSrescaled) > 0.01)
hvDS = hvPlotters.hv.Dataset(xrAFfiltered.unstack('BLM').fillna(0).squeeze(drop=True))
# hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM').fillna(0).squeeze(drop=True))
# Plot heatmap + ADMs
(hvDS.select(l=[2,4,6], m=0).to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap, clim=(-0.5,1)).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1) #.overlay(['l','t']).layout('Orb').cols(1)
In [ ]:
# Line-outs
# Note this is quite slow to render
data.BLMplot(backend='hv', xDim = 't', ylim=(-1,2), filterXS=0.01)
In [ ]:
hvDS.select(l=[2,4,6], m=0)
Versions¶
Code¶
In [ ]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [ ]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews', 'pandas'])
In [ ]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
In [ ]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
Notes¶
22/03/24
v1, quick run based on previous OCS code plus a few recent function updates.
TODO:
- Testing with alternative polarization states.
- More careful tests/checks in low-XS regions (currently just filter out, but see BLM spikes in a few cases without XS filter applied).
- More careful comparison with R-matrix results.
- AF-PAD plots, grids should work here...
Scratch¶
In [ ]:
break
In [ ]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))
#*************** Throw everything in method - breaks for large datasets (RAM issues)
# data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}) # NOTE THIS USES independently set ADMs!
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3) # OK for 61 t-points (x2 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3, thresDims = ['Eke','t']) # NEED TO ADD DIM CHECKING TO thresDims
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-4) # 1e-4 working OK for 30 t-points (x4 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-5) # 1e-5 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-6) # 1e-6 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=None) # FAILS for 16 t-points (x8 downsample) at 60Gb RAM on Jake (Swap file size issue????)
# OK for 8 t-points (x16 downsample) at ~55Gb RAM peak on Jake.
# thres=1e-4) BREAKS KERNEL #, thresDims=['Eke','t']) # TODO: currently thersDims must be in matE only for AFBLM routine.
# 21/07/22 version - above only gives L=2 max, missing L=4,6. But bumping threshold to 1e-4 kills calculation...
# data.AFBLM(keys = keys, selDims = {'Type':'L'}) # Default case... sets/uses ADMs per key
# UPDATE: testing with scale factor for ADMs, and forcing to real - might be better?
#*********** Using grouping to avoid memory issues
# QUITE slow with current code, but more robust
# TODO: make this better with basis set recycle (cf. fitting code)
#
import pickle
def calcAFBLMGroup(data,ADMs,keys,AFBLMdict):
"""Wrap for XR groupby"""
# Class routine - returns to main datastruture
data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}, thres=None) #.map(standardize)
# Copy to new structure & restack
# keys = data._keysCheck(keys)
# AFBLMdict = {}
# Restack data to lists for merge
for key in keys:
AFBLMdict[key]['AFBLMlist'].append(data.data[key]['AFBLM'].copy())
return AFBLMdict
AFBLMdict = {k:{'AFBLMlist':[]} for k in keys}
# tBins = np.arange(38,39.1,0.5) # Test range
tBins = np.arange(trange[0], trange[1] + 0.1, 0.5)
# trange
import time
for item in data.data['ADM']['ADM'].groupby_bins("t", tBins):
start = time.time()
print(f'Calculating for group: {item[0]}, {item[1].t.size} t-points.')
AFBLMdict = calcAFBLMGroup(data, item[1], keys, AFBLMdict)
end = time.time()
print(f'Delta: {end - start}s')
# Concat XR lists
for key in keys:
AFBLMdict[key]['full'] = xr.concat(AFBLMdict[key]['AFBLMlist'],"t")
AFBLMdict[key]['full'].name = f'{key}_AFBLM'
# Push back to main class
data.data[key]['AFBLM'] = AFBLMdict[key]['full'].copy()
# Save object with Pickle
with open(f'OCS_{key}_AFBLM_220722.pickle', 'wb') as handle:
print(f'Saving {key} with pickle')
pickle.dump(AFBLMdict[key]['full'], handle, protocol=pickle.HIGHEST_PROTOCOL)
# Save full class with Pickle
# NOW GIVING ERRORS ON RERUN (although worked previously), AttributeError: Can't pickle local object 'plotTypeSelector.<locals>.<lambda>'
# NOT SURE WHERE THIS IS COMING FROM AS YET
# UPDATE: IN data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails']
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
# import pickle
with open(f'OCS_AFBLM_220722.pickle', 'wb') as handle:
data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data (full class) with pickle')